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Abstract. Our paper presents an attempt to axiomatise signal process- 
ing. Our long-term goal is to formulate signal processing algorithms for 
an ideal world of exact computation and prove properties about them, 
then interpret these ideal formulations and apply them without change to 
real world discrete data. We give models of the axioms that are based on 
GAUSSian functions, that allow for exact computations and automated 
tests of signal algorithm properties. 
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1 Introduction 
1.1 Motivation 

In signal processing we consider real or complex valued functions. These func- 
tions represent signals or frequency spectra and their arguments are considered 
to be time values or frequency values, respectively. There are some fundamen- 
tal operations like pointwise multiplication "•" and convolution "*" of signals. 
The Fourier transform F converts between signals and frequency spectra. For 
precise definitions of these operations see Section [2] Additionally there are some 
essential laws, that every signal processing scientist is familiar with, such as the 
law, that the Fourier transform is an homomorphism, that maps multiplication 
to convolution: 

T(x ■ y) = Tx * Ty. 

This is an important connection and it is amazingly simple, but it is not quite 
true. 

— The first problem is, that depending on the precise definitions of the in- 
volved operations there may be a factor 2n or \/2tt to make the above 
identity correct. This is like working in the imperial system, where a gal- 
lon is not just the cubic power of a length unit but 231 cubic inch. This 
is at least cumbersome and error-prone when done manually, but it also 
complicates computer implementations. With the transcendent factors we 
have to work in fields like Q(\/27r) or Q(tt) or would have to work with 
approximations, but working with only rational numbers would of course 
be still exact and more efficient. However, if we succeed to suppress the 
transcendental factor in the above identity by the change of a definition, 
then it will show up in another place. 
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Thus our first contribution is to select a set of operations that are com- 
mon in signal processing and list identities that show their interrelation in 
principle. We use the degrees of freedom in the operation definitions for 
simplifying the laws as much as possible, that is, avoid constant factors 
and assert maximum symmetry in Section[5] Metaphorically speaking, we 
try to define something like a metric system for signal processing. 

— The second problem is, that the convolution and Fourier transform are 
not always defined. For instance the straightforward definitions of the 
Fourier transform and its inverse as given in Section [2 map from abso- 
lutely integrable functions to bounded functions, i.e. J", J 7-1 : £i(K) — > 
Coo (K) • That is, in general you cannot apply the inverse Fourier trans- 
form to the result of a Fourier transform. There are even more such dif- 
ficulties, but they can be resolved by restriction to the Schwartz space S 
of arbitrarily smooth and rapidly decaying functions. 

— The third problem is, that commonly signal processing in a computer is 
performed on discretised functions with finite precision numbers. That is, 
laws as the above one do not hold exactly. 

Our contribution concerning this problem is in the subsections of Section^ 
There we take the presented laws as axioms and develop models in terms 
of extended GAUSSian functions. All of these functions are in S. The more 
we extend those functions by factors and parameters, the more operations 
we can support. Our goal is to give definitions of functions, such that we 
only need to cope with rational parameters. This way we can represent 
a wide range of functions, we can compute exactly and we could even 
generalise these models to any algebraic field. 

1.2 Basics 

The theory of Algebra provides many axiomatically described structures, that 
generalise common mathematical objects and operations on them. E.g. groups 
abstract permutations, lattices abstract BoOLEan logic, rings abstract integer 
arithmetic, fields abstract rational arithmetic. If we are able to perform a proof 
with the axioms of a particular algebraic structure, then our proof automatically 
applies to every such structure. In this paper we would like to abstract from what 
is commonly called signal processing. 

Signal processing is the research area of construction, transformation and 
analysis of oscillation functions in one variable (the time), approximation of 
those functions and related operations with discretised data and efficient im- 
plementations on digital computers. An important view on an audio signal is 
the frequency spectrum, because that is closer to the way humans hear than 
the time-domain representation of a signal. We obtain a frequency spectrum by 
the Fourier transform of a signal. Closely related to signal processing is im- 
age processing, with the main difference being, that functions in two variables 
are studied. Another related area are random distributions in stochastics. The 
convolution, a common operation in signal processing, of two random distribu- 
tions yields the distribution of the sum of two random variables. The Fourier 
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transform yields the characteristic function of a random distribution and there 
are more connections (see Section [3]). 

However, signal processing, image processing and random distributions are 
relatively seldom viewed from an algebraic point of view. This may have to do 
with the tradition to focus on signal values (e.g. f(t)) rather than on larger 
objects like signals (e.g. /). E.g. it is common to express the delay of a signal 
by fit — d) rather than to use an operator like in / — > d. On the one hand this 
has the advantage, that some properties are obvious (e.g. fit — [d + di)) = 
f{{t — do) — di)), since they are a consequence of simple arithmetic, but are 
not obvious at the higher level (i.e. / — > [d + di) = (f — > di) — > d ). On the 
other hand we cannot cleanly express identities involving intrinsically functional 
transforms like the Fourier transform, that has lead to custom notations like 
F[u>) <h» fit), that could be expressed the functional way as FF = f . 

2 Finding a set of fundamental operations 

We like to start this section listing the operations, that are commonly used in 
signal processing, together with typical applications. We begin with the opera- 
tions with indisputable definitions and continue with the ones, where differing 
definitions are around in the literature. For the variant definitions we show, what 
laws they imply. Then we choose the definitions that make the laws most sim- 
ple. We close the section with a comprehensive list of laws that hold for our 
definitions. 

In Table [T] we list the signal operations together with their definitions and in 
Figure [1] we illustrate them using example signals. 

— Shrinking a signal as in ([lj means to increase all contained frequencies 
proportionally and to shorten time accordingly. 

— Translating a signal as in @ means to delay it. 

— Summing two signals as in ([3]) means to superpose them, that is, to play 
them simultaneously. 

— Multiplying two signals as in (j4]) means to control the amplitude of one 
signal by the other one. For certain choices of signals this is also known 
as ring-modulation. 

— Convolving two signals as in ([5]) means to apply the sound characteristics of 
one signal to the other one. It may be used for suppressing or emphasising 
certain frequencies, for smoothing or for application of reverb. 

We continue with the controversial definitions. Actually, there is essentially 
one operation, that is defined differently throughout the signal processing liter- 
ature: The Fourier transform. It is the transform that computes the frequency 
spectrum for a signal in the time domain [backward or analysis transform) and 
vice versa [forward or synthesis transform). For an example see ([6]) in Figure [TJ 
The other definition with varying instances is the one for functions that are 
eigenfunctions of the Fourier transform. Obviously it depends entirely on the 
definition of the Fourier transform. 
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Table 1. Basic signal processing operations 



operation 



application 



definition 



shrinking 

translate delay 
adjoint 

sum mixing 

multiplication envelope 

convolution frequency filter 

modulation frequency shift 



alter pitch and time (x l k)(t) — x(k ■ t) 



(x ->■ k)(t) = x(t - k) 

x* — x J, — 1 
(x + y)(t) = x(t) + y(t) 
(x ■ y)(t) = x(t) ■ y(t) 

(x*y)(t)= / x{t) ■ y(t - t)At 

x ■ cisl where cisl t — exp(27ri • t) 



12 = 




(1) 



(2) 



(3) 



(4) 



(5) 



(6) 



Fig. 1. Illustration of the basic signal processing operations. 
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1. Oscillations with period 1 

Fourier forward T\x{t) = I cxp(27rz ■ r • t) ■ x(t) dt 

Jr 

Fourier backward J"f 1 x(t) = / cxp(— 2iri ■ r • t) ■ x(t)dt 

Jr 

duality Fi{F\x) = x \. — 1 

eigenfunction g(t) = exp (— tt ■ i 2 ) (7) 

eigenvalue is 1 T\g = g 

unitarity {x, y) = (Tix, Fiy) 

convolution theorem F\{x * y) — T\x ■ T\y 
T\{x ■ y) = Fix * T\y 
derivative T\{x'} = — 2m ■ id -T\x (8) 

2. Oscillations with period 27T 

Fourier forward Tix(t) = — ^= ■ j exp(i • r • t) ■ x(t)dt 

v27r Jr 

Fourier backward T^ 1 x(t) = —== ■ / cxp(— i ■ r • t) ■ x(t)dt 

2 ^ Jr 

duality Ti{Tix) = x I — 1 

eigenfunction g{t) = exp (— t 2 ) (9) 

eigenvalue is 1 Tig = g 

unitarity (x, y) = {T^x.Tiy) 

convolution theorem Tiix * y) = V^n ■ TiX ■ Tiy (10) 
^(x-y) = — ^= ■T 2 x*T- 1 y (11) 

V 27T 



derivative Ti(x'} = —i ■ id-J^x (12) 

3. Oscillations with period 27r and no roots of ir 

Fourier forward Fsx{t) = / exp(z • r • t) ■ x(t) dt 

Jr 

Fourier backward J 7 3 ~ 1 x(t) = — ■ / cxp(-z ■ t ■ t) ■ x(t) dt 

2?r J r 



We do not further consider Definition 3, because forward and backward 
Fourier transform have different factors, thus laws for forward and backward 
transform differ in factors. For definitions 1 and 2 many laws are equal up to the 
choice of the transform direction. 

At the first glance the definitions 1 and 2 of the Fourier transform seem 
to be equally convenient or equally inconvenient. Thus many textbooks just in- 
troduce a definition and do not explain, why they prefer the one to the other 
possible ones. However, we think that factors in laws are bad, especially worse 



6 Henning Thielemann 



than factors that can be hidden in a function definition. E.g. compare the eigen- 
function property of the GAUSSian function in Q and ©: We can define the 
GAUSSian bell curve in both ways and call it just g. The factor ir does no longer 
get in the way when transforming equations containing T and g. The same ap- 
plies to the laws on derivatives in (J5J and (TT^I) , where we can consider 2ni ■ id as 
one function. In contrast to that, the convolution theorems in (TTU)) and (fTTj) have 
constant factors, even different ones. This makes manual equation manipulation 
cumbersome and error-prone. Even more it makes computations exclusively with 
rational numbers impossible, due to the irrational factor \p2/K. We could suppress 
the factors in (|10[) and by defining convolution or multiplication containing 
a constant factor, but we think this is not natural. 

For these reasons we will stick to definition I and omit the index of T in 
the rest of this paper. Below we list the laws that follow from these definitions 
and that we want as axioms for our algebraic structure. Be aware, that we have 
omitted convergence conditions for the operations that involve integration. The 
laws marked with # can be derived from the remaining laws and are just given 
for the purpose of completeness. 



x + y = y + x 
x ■ y = y ■ x 

x*y = y*x^ 



x^O 
x I I 



-> a = 



(x + y) 
(x-y) a 
(x*y) a 

(x ■ y)' 



X 
X 



(x — > a) + (y — > a) 
(x a) ■ (y a) 
x * (y a) 
x • y + x ■ y 



x + (y + z 
x- (y z 
x ■ (y + z 
x * (y * z 
x * (y + z 
(x — > a) — > b 
(x I a) I b 
(x lb) — > a 
(x + y) I a 
(x-y) la 
(x * y) I a 
(x * y)' 



(x + y) + z 

(x-y) ■ z 

x ■ y + x ■ z 

(x * y) * z ^ 

x*y+x*z^ 

x — > (a + b) 

x I (a ■ b) 

(x -+(a-b))\rb 

(x I a) + (y I a) 

(x I a) ■ (y I a) 

\a\ ■ (x I a) * (y I a) 

x * y' 



T(x + y) = Tx + Ty 
T(x * y) = Tx ■ Ty 
(x, y) = (Tx, Ty) 

T(x la) = \a\ ■ Tx I — 
a 

T(Tx) = x I — 1 
T(x*) =Tx 



T(k ■ x) = k ■ Tx 
T(x ■ y) = Tx * Ty & 

IMIa = W^h * 
T(x — > a) = (cis I a) ■ Tx 



T(x') = —2iri ■ id -Tx 
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3 Development of Gaussian models 

Now that we have stated some axioms, we want to construct an ideal world, 
that is, a class of functions (or signals), where they hold. This class of function 
shall allow for simple constraints of the laws, for exact and efficient computa- 
tion (more precisely: operations in a field), and shall allow to represent many 
mathematically important objects exactly and real world signals approximately. 

Since GAUSSian functions are eigenfunctions of the Fourier transform, they 
are perfect for representing signals in both time and frequency domain. They can 
be extended in a relatively simple way, such that all of the operations can be 
performed, that we listed initially. 

We like to stress that GAUSSian functions are not only interesting, because 
they allow to perform the operations we want, but GAUSSian functions are cen- 
tral to signal processing and stochastics. 

— GAUSSian functions are used as filter window for smoothing. 

— GAUSSian functions "minimise uncertainty" , that is, they are the functions 
where Heisenberg's uncertainty relation becomes an equation. [7] 

— With complex modulation GAUSSian functions are called Gabor atoms 
or time-frequency atoms in the Gabor transform. The Gabor transform 
is a windowed Fourier transform, that shows how frequency components 
evolve over time. 

— Complex modulated GAUSSian functions with a correction offset are called 
Morlet wavelets and are used in the Continuous Wavelet Transform. 
This transform is also intended for showing the evolution of frequency 
components over time, but it has higher time resolution and less frequency 
resolution for high frequencies. 

— Best basis pursuits and matching pursuits are techniques for decomposing 
a signal into a finite number of irregularly located time-frequency atoms. 


Since Gabor transform, Morlet wavelet transform as well as best basis and 
matching pursuits aim at decomposition of a signal into time-frequency atoms, 
we have several tools for approximating real world signals using those atoms as 
building blocks. 

In stochastics the density of the Normal distribution is a GAUSSian function. 
The Central Limit Theorem states, that adding more and more random variables 
(and divide by the square root of the number of added variables) in most practical 
cases approaches a normally distributed random variable. Translated to signal 
processing this means, that smoothing a signal again and again with practically 
any filter window, approximates a GAUSSian filter. 

The derivation of the representations below is not particularly difficult, but 
we focus on finding representations that simplify most operations in terms of 
use of transcendental constants and irrational algebraic functions. We have im- 
plemented these function classes in Haskell using the NumericPrelude type class 
hierarchy. You find our implementation at 

http : / / code .haskell . org/numeric-prelude/ src/MathObj /Gaussian/ 
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3.1 Simple Gaussians 

Simple GAUSSian functions shall be the first class of functions that we want to 
consider. In order to avoid a transcendental factor containing tt in any of our 
function parameters when applying Fourier transform J 7 , we cannot choose 
f(t) = cxp(— t 2 ) but we have to choose its eigenfunction f(t) = cxp(— 7r • t 2 ). We 
also want to support shrinking of a function, what requires adding a shrinking 
parameter c, that we like to write in curried form: f(c)(t) = exp(— tt ■ (c • t) 2 ). 
However this would yield a square root in the convolution of two such functions. 
It can be prevented by choosing the form f(c)(t) — cxp(— tt-c4 2 ) with c € Q. The 
Fourier transform of a shrunk function leads to another factor, the amplitude 
y of the function. We end up with the function form: 

f(y, c )(t) = Vv ■ cxp(-7r • c • t 2 ) 
ceQ y e [0,oo) n Q. 

This simple function class already allows for several operations, where con- 
volution and Fourier transform have constraints that assert convergence: 

scaling 
shrinking 
conjugate 
multiplication 
power with r > 

convolution c + c\ > 

Fourier transform c > J 7 ^ 1 (f(y, c)) = / 

Typical functional like function norms do not easily satisfy our goal of using 
the most simple algebraic structures. They need roots of parameters or constant 
transcendental factors, and thus require special treatment: 

£i-norm 



£ 2 -norm 

£oo-norm 

£„-norm 



k ■ f(y,c) 
f(y,c) I k 

7M 
/(yo,c ) • /(yi,ci) 
f(y,c) r 

/(yo,c ) * /(yi,ci) 



f(yk z ,c) 
f{y,c-k 2 ) 
f(y,c) 

f{yo-yi,co + ci) 
f(y r ,r- c) 

j. ( yo ■ yi c ■ ci \ 

V C + Ci ' C + Ci / 




\\f(y,c)\\ 2 
H/(i/,c)lL 
\\f(y,c)\\ p 

\t^t 2 ■ f(y,c)(t)L 



'2-c 
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3.2 Translated and modulated Gaussians 

In the next step we want to translate and modulate the GAUSSian function. The 
most simple function class, that allows this, seems to be: 

f(y, a, b, c)(t) = sjy • cxp (-7T • (a + b ■ t + c ■ t 2 )) (13) 
y e [0, oo) n Q 
ceQ {a,b}cQ + iQ. 

With this representation we can perform the following operations: 

translation f(y, a, b, c) — > k = f(y, a — b ■ k + c ■ k 2 ,b — 2 • c • k, c) 

modulation f(y, a, b, c) ■ (cisl \.k) = f (y,a,b + 2 • i ■ k,c) 

scaling k ■ f(y, a, b, c) = f(y ■ k 2 , a, b, c) 

shrinking f(y, a,b,c) \,k = f(y, a,b ■ k,c- k 2 ) 

conjugate f(y, a, b, c) = f(y, a, b, c) 

power with n > f(y, a, 6, c) n — f(y n , n ■ a 7 n ■ b,n ■ c) 

multiplication 

f{yo,a ,bo,co) ■ f(yi,ai,bi,a) = f(y ■ yi, a + oi, b Q + bi, c + ci) 

convolution c + c\ > 

f(yo,ao,bo,c ) * f(yi,ai,bi,ci) = 

f ( Va-yi „ n i _ _ (fco-frl) 2 bo-cx+bx-cg c -ci A 
J ^co+ci'" OTU1 4-(c +Ci) ' Co+ Cl ' Co+Ci J 

Fourier transform c > 



T-\f{y,aAc)) = f^a- b -,- 1 -^^ 



The correctness of the equations for the Fourier transform and the convolution 
are not so obvious. The Fourier transform can be derived by translating the 
function to the origin and demodulate it, such that it becomes real. Then do 
Fourier transform and translate and modulate it corresponding to the normal- 
isations that we performed in time domain. The convolution can also be derived 
from such an normalisation. An alternative is to multiply in frequency domain. 
We could also employ the definition 

f(y, a, b, c)(t) = ' exp(-(a + b-yfirt + c- Trt 2 )) (14) 

and the formulas for most operations would remain the same, but translation, 
shrinking and modulation would have to be interpreted with respect to the 
unit \fii. 

With the considered representation we can also represent two other kinds of 
functions that are important to signal processing: A decaying exponential curve 
can be obtained with c = A b > 0. It is frequently encountered as envelope of 
percussive sounds. Unfortunately that curve is unrestricted in time, whereas in 
natural sounds the envelope usually starts suddenly somewhere in time. 
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The other important signal is a tone of linearly changing frequency, a chirp. 
In order to represent it, we have to generalise the real parameter c to a complex 
parameter. Setting 3(c) ^ 0, 5ft(c) = yields a pure chirp, whereas 5R(c) > 
yields a chirplet, that is a chirp, that is localised in time. The name is chosen 
analogously to wavelets. The Fourier transforms and norm computations con- 
verge, if and only if, Sftc > (i.e. only for chirplets), and the convolution exists 
if and only if 5R(co + ci) > 0. The Bluestein transform allows us to write a 
Fourier transform in terms of a convolution with a chirp. 

T- X x = ((x ■ /(l, 0, 0, -»)) * /(l, 0, 0, *)) • /(1, 0, 0, -i) 

The downside of this generalisation is, that we need to maintain a complex 
amplification factor y. This involves a complex square root and we must choose 
the right branch. For a single convolution or FOURIER transform choosing the 
branch with positive real part is just the right thing. The real part cannot vanish, 
since then the integrals do not converge at all. But when multiplying square roots 
with the convention of positive real parts then we must respect 

\/co • \[c\ = ' V c o ■ ci 

!1 : (upper Co A upper ci A ^upper(co • ci))V 
(-i upper Co A -i upper c\ A upper(co • ci)) 
: otherwise 
upper c = 3c> V (3c = A 5Rc < 0). 

That is we must maintain the sign of the real part separately. In the representa- 
tion (|13[) we can implement a flipped sign by adding i to b. However maintaining 
this sign means comparisons and thus would not work for generalisations from 
Q to other fields, e.g. finite fields. 



3.3 Gaussians multiplied with polynomials 



Another important operation in signal processing is derivation, be it in differ- 
ential equations like oscillation equations, for the representation of an eigen- 
basis of the Fourier transform, or as a highpass filter, that is, a frequency 
filter that emphasises high frequencies and suppresses low frequencies. Deriva- 
tion requires, that we extend our representation to a product of a GAUSSian 
function and a polynomial function. However using the simple representation 
f(y, a, b, c)(t) = ,/y ■ exp (— tt • (a + b ■ t + c ■ i 2 )) from (TI3)) , this would mean to 
maintain polynomial expressions of n as coefficients of the polynomial factor. 
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We like to avoid that and instead extend (fT^Ij) to: 

f(y, a, b, c)(t) = y/y ■ cxp (-(a + by/n ■ t + c ■ ir ■ t 2 )) 
tp((y,a,b,c),p)(t) = f(y,a,b,c)(t) ■ p(t) 
ye [0,oo)nQ 

p(t) G (Q + iQ)[V^-t] 

n 
3=0 

As mentioned above translation and modulation have to be interpreted with 
respect to a unit y/w, and derivation must contain a factor of ypix . However in 
order to get the usual units we can still replace Q by Q(y/n). 

For the following list of instantiations of signal processing transforms we like 
to subsume the parameters of / in a parameter tuple a. 

translation ip(a,p) — > y/irk = (fa — > y/nk) ■ (p — > y/nk) 

modulation ip(a,p) ■ (cisl 4- —p=) — ( foi ■ (cisl I —= ) ] -p 

scaling fc • <p(a,p) — fa ■ (k ■ p) 

shrinking ip(a,p) I k = (fa I k) ■ (p ! k) 

conjugate <p(a,p) — fa-p 

multiplication ip(a Q ,p Q ) ■ cp(a 1 ,p 1 ) = (f(a Q ) ■ f(a%)) ■ p^P\ 
power with neN ip(a,p) n = (fa) n ■ jf L 

convolution (p(a ,p ) * cp(a 1 ,px) = T (J" _1 (^(q;o,Po)) ■ {vip-uVx)^) 

Fourier transform 

F-\<p{a, s:p)) = s- F~ l {fa) + ^= • (F-\y(a,p)))' 

where sTp(t) = s + y/ir ■ t ■ p(t) 

Differentiation 

-j= ■ (cp((y, a, b, c),p))' = f(y, a, b, c) • {-= ■ j? - (t ^ b + c ■ yfH ■ t) ■ pj (15) 
Integration 

T 

■ Loo p(0/> a > 6 ' °)^p)( t ) dt 



s ■ exp —a 



f(y,a,b,c)(t)dt + f(y,a,b,c)(T)-q(T) 
b 2 \ 1 + ^(tTc + ^- t 



Ac J 2yfc 

+f(y,a,b,c)(T)-q(T) 

l.q>( t )-p(t)-s 
where q(t) = — — — (16) 

yW b+C-y/^-t V ' 
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Eigenfunction of Fourier transform 

e n = /(l,0,0,2)(">-/(l, 0,0,-1) (17) 

The equation (fllj)) is the inverse of (IT51) . This implies that in (TIB)) the polynomial q 
depends recursively on itself. However because the degree of p is one more than 
that of q, the leading term of q only depends on the leading term of p. Thus we can 
successively determine the terms of q starting at the highest one. The equation 
can be translated almost literally to a polynomial division with remainder s in 
our Haskell implementation and be solved using lazy evaluation. 

In (IT7|) we have used the definition of Hermite polynomials and the known 
fact, that the GAUSSian function multiplied with a Hermite polynomial is an 
eigenfunction of the Fourier transform. 

3.4 Mixed Gaussians 

In order to support sums of signals, we must maintain a set A of parameters for 
the GAUSSians and a map P from the GAUSSian parameters a to the associated 
polynomial factor. 

fa ■ Pa 

Eventually, this representation is general enough in order to be target of a win- 
dowed Fourier transform or a best basis or matching pursuit. That is we can 
approximate real world signals in a natural way and perform exact signal pro- 
cessing operations on them. 

4 Related work 

With our paper we wanted to draw a connection between Computer Algebra on 
the one side and Signal Processing and Stochastics on the other side. With "Com- 
puter Algebra" we mean exact computations involving complex mathematical 
objects like polynomials, polynomial ideals, groups, that is at a higher level than 
computing with individual numbers but at a lower level than computing with 
general mathematical expressions as in symbolic manipulations. Although signal 
processing is certainly not the most prominent application of computer algebra, 
there are many problems that were solved using computer algebra methods. [4 
The most famous application is probably the Discrete Fourier Transform, that 
can be performed in log-linear time with respect to the length of the input data 
using techniques from number theory, finite fields, polynomial rings and auto- 
mated code generation. |1I3| Closely related is the fast convolution of discrete 
signals, that uses the Fast Fourier Transform and by the chirp transform it is 
also possible to express a fast Fourier Transform in terms of fast convolution 
algorithms. Another computer algebra application in signal processing is the 
design of frequency filters, where we have to construct rational functions given 
conditions for the location of its zeros and poles. [5] 
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In [5] the authors develop signal processing the algebraic way, as we do in 
our paper. That is, opposed to the sample value focus of most signal processing 
literature they treat signals as objects, define operations on them and propose 
and prove laws. The book is concerned with two-dimensional signals, but the 
difference to one-dimensional signals is not essential in this approach. Unlike our 
consideration of real functions and especially modified GAUSSian functions the 
authors stick to discrete signals. 

Compared to established computer algebra systems and their symbolic inte- 
gration machineries, our framework provides no new class of closedly integrable 
functions. All of Maple 10, Mathematica 5.2, Maxima 5.20.1 can integrate the 
integrals occurring in convolution, Fourier transform, norm and scalar product 
of our products of Gaussian function and polynomial in closed form. MuPAD 
4.0.6 cannot integrate the more complicated convolutions and Axiom 20091101 
cannot cope with the integrals at all, since it does not yet support assumptions. 
We need assertion for the coefficient c of the quadratic term in the exponent of 
the GAUSSian function. It must have positive real part, or must be at least a 
positive real number. That is, with an assertion we must exclude signals with 
constant amplitude or even unbounded amplitude. 

5 Outlook 

5.1 Dirac impulses 

The Dirac impulse S is a virtual function that is infinitely high at the origin 
and zero elsewhere, enclosing an infinitely high and infinitely narrow rectangle 
of area 1. If this function would exist, it would be the identity element of con- 
volution. Its Fourier transform is the function that is constant 1, because this 
is the identity element of pointwise function multiplication. In stochastics the 
Dirac impulse is needed for representing mixed discrete/ continuous probability 
distributions. Existence of derivatives of the Dirac impulse would eliminate the 
need for a distinct differentiation operation, because x' = x * 8'. We could also 
more easily hide the factor ^/tt in the differentiation operation and we could 
represent frequency spectra of polynomial functions. 

Several approaches like Schwartz distributions and non-standard analysis 
were developed, in order to get a strictly founded notion of a Dirac impulse. 
However, none of them is completely satisfying: Schwartz distributions have 
no longer a notion of function application and they cannot be multiplied point- 
wise. Non-standard analysis allows to define infinitely high and infinitely narrow 
functions, that actually let real functions unaltered at the "coarse" scale of real 
numbers after convolution. However when considering a non-standard function 
on all scales, convolution with the non-standard Dirac impulse well changes the 
convolution partner. It is not known to us, whether an approach can exist at all, 
that fulfils all expectations. 

All the more it is interesting whether we can have an object, that exactly 
behaves like a Dirac impulse in our theory. Since our theory is abstracted from, 
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but not bound to real functions, we could check this way, whether a DlRAC 
impulse makes sense at all. Formally in our approach a DlRAC impulse could 
be represented by /(l, 0, 0, +oo). The term +oo could be made precise by using 
projective geometry, i.e. by allowing an object like q. But then it is open, whether 
we should use 

1. ex P (-7T- »+b-t+c-t 2 y 

2. cxp (-7r ■ (a + b- 1 + | • t 2 )), 

3. exp ^-7T • (a + b ■ | + c • ^fj or 

4. ex p(-^ + ^.t+^.t 2 ) 

with a projective interpretation of the fractions, and how to cope with the am- 
plitude parameter. 



5.2 Discrete signal processing 

We would like to have the same set of operations and laws for discrete signals 
that we already have for real signals. In order to have dual time and frequency 
domains, we need to content ourselves with periodic discrete signals. For instance 
for discrete periodic signals x and y of period length n convolution and Discrete 
Fourier Transform (DFT) are usually defined as: 



{x*y)k = 



Xj ■ yk-j 



( DFT ^ = ^-£ ex p(~;-fc 

(DFT -1 x)k = -4= • exp (—^- ■ j ■ k 
v" hi V n 



j£2 

The factor is chosen, such that DFT becomes unitary. However with this 
definition it does not hold DFT(x-y) = DFT a; * DFT y, but instead ^/E-BFT(x- 
y) = DFT x*DFT y. One solution would be to add the factor to the definition 
of the convolution. This is at least very uncommon. An alternative is to treat 
discrete signals as piecewise constant functions. The sums are turned to integrals 
and thus need a step width. To this end we equip every signal with a sampling rate 
and denote it with rate. It holds rate(DFTa;) = " . We obtain the definitions 

^ ' rate x 



{x * y)k = — - — • ^ Xj ■ yk-j for {rate x, rate y} = {rate(x * y)} 

jez„ 

(DFT^.-l-.^exp^.^fc).^ 
(DFT-,), = ^-E «P (~- 
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Following the reasoning of real signals, we need eigenfunctions of the Fourier 
transform. Generic simply representable eigenbases of the Discrete Fourier 
transform are currently not known, but according to the PoiSSON summation 
formula, the discretised and periodically summed eigenfunctions of the Contin- 
uous Fourier transform, are eigenvectors of the discrete transform. However 
discretising an eigenbasis of the continuous Fourier transform may not yield a 
discrete eigenbasis. 

In discrete signal processing the identity element of convolution is simple 
to get: It is the signal that is 1 at index and zero elsewhere. In contrast to 
that, the operation of signal dilation may lead to undefined and multiple times 
defined elements in the resulting vector. The natural solution is to set undefined 
elements to zero and sum up all candidates for multiply defined output elements. 
This can be written generally in the following way, where n is the signal period 
length and the empty sum is zero: 

0t%= Xl ~ 

Z6Z„: l-k=j 

This definition matches shrinking the vector in the frequency domain. Nonethe- 
less, we have to drop invertibility of dilation from the list of laws, that hold for 
real signals. 

Another problem is the definition of differentiation. We could replace it by 
centred discrete differences. Via the Fourier transform this would also yield a 
notion of periodic polynomials, namely polynomials in sin instead of t. But 
this interpretation of differentiation is different from differentiation of a con- 
tinuous function with subsequent discretisation and periodic summation. Thus 
it cannot be used for eigenvector computation in the same way we used it for 
continuous signals. 

Summarised, we cannot simply perform the operations on our parameter 
tuples, that we developed for continuous signals, and use them for discrete signals 
by just interpreting them in a discrete way. 
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